Potential for regulatory genetic networks of gene expression near a stable point 
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A description for regulatory genetic network based on generalized potential energy is constructed. 
The potential energy is derived from the steady state solution of linearized Fokker-Plank equation, 
f — , and the result is shown to be equivalent to the system of coupled oscillators. The correspondence 

between the quantities from the mechanical picture and the steady-state fluctuations is established. 
Explicit calculation is given for auto-regulatory networks in which, the force constant associated 
\ with the degree of protein is very weak. Negative feedback not only suppresses the fluctuations but 

! * , also increases the steepness of the potential. The results for the fluctuations agree completely with 

O ' those obtained from linear noise Fokker-Planck equation. 

O: 

A regulatory network of gene expressions consists of a group of genes which co-regulate one another's expressions. 
Such networks provide a fundamental description of cellular function at the DNA level. Recently, the advance of 
experimental techniques in constructing synthetic networks with the ability of monitoring them has provided some 
essential elements, such as switchp], @[ Q and oscillator [i], for the design of biological circuits. In modeling 
|| the dynamics of a regulatory network, rate-equation approach is often used; the approach reflects the macroscopic 
JL ■ observation with deterministic nature. However for systems with small molecular number, intrinsic fluctuations 
. ' become important. The noise- induced effect may be incorporated into the framework by employing the master 
, equation and then proceeding via stochastic Monte Carlo simulations. In general, master equation is discrete in 
nature. By using the technique of O-expansionQ, we may convert master equation to continuous Fokker-Planck 
equation which, then, is managed analytically by various approximations. Significant progress has been made along 
this line in understanding the regulation mechanism. One of the noticeable examples is the auto-regulatory networks 
of a single gene for which, the protein, encoded in the gene, serves as the regulator of itself through either negative 
or positive feedback. Such autoregulation is a ubiquitous motif in biochemical pathways. It was demonstrated by 
Becskei and Serrano that an autoregulatory network with negative feedback may gain stability Q ■ Further analyses was 
given by Thattai and van OudenaardenQ and by Ozbudak et al.@, and the results indicate that noise is essentially 
determined at the translational level and negative feedback can suppress the intrinsic noise. Moreover, Tao and Tao 
et al. used the linear noise Fokker-Planck equation to study the fluctuations and obtained the results consistent 
qualitatively with previous works [Tol HlT|. 

One may conclude from the results above that the intrinsic noise associated with a genetic network is closely related 
. to its regulation scheme. This Letter then attempts to provide a physical picture on this relation via the establishment 
of a mechanical analogous system. To achieve this, we first construct the solution of non-equilibrium steady state 
for the Fokker-Planck equation near a stable point. Then, the potential of the system, defined as the negative of the 
1 logarithm of the solution, can be first approximated as an harmonic oscillator potential . Subsequently, we introduce a 
measure for the steepness of the potential near a stable point and give the exact relations between the force constants 
of coupled oscillators and the correlations of fluctuations. Thus, the physical property of a regulation scheme can 
be revealed from the corresponding mechanical analogue specified by the force constants of coupled oscillators. This 
paper starts with the general construction for a ei-dimcnsional regulation network, followed by the explicit calculations 
^ ' of auto-regulatory networks. 

Consider a <i-dimensional regulatory network of gene expression. The network is specified by the macroscopic rate 
equations 
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Xl = fi 0) (1) 

with i = 1, 2, d, and the drift force fi defined as 

f i (x) = R i (x)-fax i . (2) 

Here, x T = {x\, X2, Xj) with the superscript r for the transpose of a vector, Xi represents the concentration of 
mRNA or protein, the function Ri (x) describes the synthesis or feedback regulation of molecule i, and the constant 4>i 
denotes the degradation rates of Xi. The network is assumed to form a chain with the nearest neighboring regulation, 
Ri (x) = Ri (xi—i, Xi + i); however, the formulation presented in this work can be extended to more complicated cases 
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straightforwardly. The fluctuation may be incorporated into Eq. JT]) by means of the master equation approach. For 
this, we introduce the volume factor Q to give the molecular number n T = n 2 , n<j) as n T — Qx T . In terms of 
molecular numbers n, the corresponding master equation of Eq. |T]) can be written as 



OP (n, t) 
Ft 



- 1) [(0<ni) P (n, t)] + Ri (*) [Ei- - 1] P (n, t) : 



i=l i=l 

where P (n, i) is the probability distribution, and the step operators Ei± are defined as 

E i± G( ni ) =G(ni±l) 



(3) 



(4) 



for a function of molecular numbers G (n). Then, the technique of Sl-expansionQ is employed to transfer the discrete 
process of Eq. ^ to a continuous process described by the Fokker-Planck equation, 



dp (x, t) 
dt 



+ V • J(x,t) = 0, 



(5) 



where (V) r = (d/dx±, d/dx2, d/dxd), p(x,t) is the distribution density, J (x,t) is the density current defined as 

1 



J {x,t) = f(x)p(x,t) 
and the elements of the diffusion matrix D (x) are 

Dij (x) = Sij 



n 



[D (x)-d]p(x,t) 



Ri (x) 



(6) 



(7) 



with the Kronecker delta Sij — 1 for i — j otherwise 0, note we do not sum over repeated indices. 

We are interested in the behavior of p (x, t) for the region near a equilibrium stable point of Eq. ((T|), say x* . After 
expanding the density current J (x,t) of Eq. © around the stable point, we obtain the linearized Fokker-Planck 
equation for the new variable y = x — x* as 



dp L (y, t) 

dt 



+ V- J L (y,t) = 0, 



where Jl (y, t), which contains only the leading order of J (x, t), is 



J L (y,t) 



F(x*)-y~^D(x*)-W 



pl (y,t) 



(8) 
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with the force matrix F (x*) defined as Fy- (x*) = dfi (x) /dxj\ and noting the fact that y itself is of order ~. 
This leads to a Ornstein-Uhlenbeck process in which, the drift force is linear and the diffusion is given by a constant 
matrix [ri IT^|. The stationary solution of Eq. ([8|), characterized by the condition V • Jf (y) — 0, can be expressed as 



with 



and 



Pl (V) = 7? exp [-$ (y)} 



/oo poo / d \ 

•••/ \l[dy m exp[-*(y)] 



*(y) = -^y T -u(x*)-y, 



(10) 



(11) 



(12) 



where Z can be referred as the partition function, and U is a real symmetric d x d matrix Note that the 
temperature in this work is always set to be 1, fc^T — 1, and hereafter we drop the arguments for all matrix elements 
known to be functions of the equilibrium stable point x* . One can determine the matrix U by substituting Eq. (|10[) 
directly into the condition V • j£ (y) = to obtain 



tr(F + -D ■ U) - y T ■ (U ■ F + -U ■ D ■ U) ■ y = 0. 



(13) 
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To solve this for the matrix U, we follow an elegant method proposed by Ao[l3| and Kwon et al.[l4| to factorize the 
force matrix as 



F 



[D + Q]- U, 



(14) 



where D is the symmetric diffusion matrix given by Eq. ([7]), and Q is an antisymmetric matrix which has to be 
determined. Such factorization amounts to decomposing the density current into two parts, j£ (y) — j? (y) + j~! (y). 
The first term of Eq. (|14[) corresponds to the dissipative part which generates a motion towards the origin with 
vanishing density current, (y) = 0; meanwhile, the second term is the cyclic part with a divergence-free current 
density, jf (y) = — (l/f2)Q • U ■ yp s L (y), which generates a circulating motion around the constant surface of $ (y). 
By substituting Eq. (jT4j) into Eq. (fl3]) . we obtain the relation 



F ■ Q + Q ■ F T = F ■ D — D ■ F T , 



(15) 



which gives enough conditions to determine the matrix Q completely. Thus, the function $ (y) of Eq. (|12p , which 
can be referred as the potential energy for the system near a stable point x* , becomes 



n 



(D + Qy 1 -F 



(16) 



A more intuitive physical picture about the characteristics of the system may be revealed by mapping <£> (y) to the 
potential energy of the system of coupled oscillators, 



*(l/)=2 



*iVi + K ij (Vi -Vjf 



i>j 



which can be castcd in the form of 



(17) 



(18) 



Note that though the regulations of the network only come from the nearest neighbors, the couplings of oscillators 
may not be restricted to the nearest neighbors. The force constants, K and k c , can be specified by equating Eq. I|18p 
to Eq. (|16|) . and the characteristics of the network near a stable point can be expressed in terms of the force constants. 
Firstly, based on the partition function of Eq. (TTTj) , which is reduced to 



Z 



d/2 



[detV (k,k c )P 1/2 



(19) 



with det V (k, k c ) for the determinant of the matrix V (k, k c ) of Eq. (|18p. we may introduce the effective free energy 
difference, AG = —hiZ, to describe qualitatively the steepness of the potential. A stable point with larger AG 
value is more easy to focus with less fluctuations. Furthermore, the variances and covariance of x\ and X2, defined as 
a i,j — ( x i x j) ~ x *i x *j> can be evaluated by using the distribution pf (y), 



Oi A = — 



: [ dy n 



ViVi exp 



-^y T -V(k, k c ) -y 



(20) 



The formulation is applied to two-dimensional regulatory networks, and the results are given explicitly in the follow- 
ings. 

Consider the case of d = 2 with regulation functions R\ (X2) and R2 {x\). The force matrix is 



F 



(21) 



where r± and ri are defined as t\ — dR\ (x^) I d x 2\ X2=x * and r-i = 9i?2 ( x i) /d x i\ Xl =x*- Then, the antisymmetric 
matrix Q, determined by Eq. (|15[) . is 



Q = 



w 
-w 



(22) 
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with W — [r^Dii — riD 22 ] / (<fii + 4> 2 )- The F and Q matrices given above with the D matrix of Eq. ||7J) yield the 
potential energy of Eq. (TTTj) as 

= § [«u/i + K 2y 2 2 + k; j2 (2/1 - y 2 ) 2 ] , (23) 

where the force constants are K\ = [—r 2 Dn + (2(f> l - r\) D22 + (2r 2 — 4> 2 + <M W] /2 (DnD 22 + W 2 ) , 
k 2 = [(2 ( /. 2 -r 2 )D 11 -r 1J D 22 -(2r 1 + 9 !. 2 -0 1 )^]/2(Li 11 D 22 + T^ 2 ), and k£ 2 

[^1^22 + r 2 Du + (0 2 — X ) W] /2 (DuD 22 + FY 2 ). For the effective free energy difference, we rescale AG by 
adding a volume factor, AG = —hiZ — ln[f2/ (27r)]; then, AG becomes half the logarithm of detV (k, k c ), and 
it is AG = (i) In [k\k 2 + (ki + k 2 ) k\ 2 ] . Moreover, for the potential energy of Eq. (f2"3")l the variances and 
covariancc of Eq. (20]) become ct 2 ! = [(k 2 + kJ 2 ) /fi] exp (-2AG), cr§ 2 = [(«! + Kf 2 ) /Q] exp (-2AG), and 
o\ 2 — \k\ 2 /f2) exp (— 2AG). Thus, expressed in terms of mechanical quantities we summarize the features of the 
system implied by the potential energy as follows. In general, the AG value characterizes the global steepness of 
the quadratic potential; the increase of the AG value makes the potential more sharper and, hence, reduces the 
fluctuations. However, for the same AG value the details of potential shape may have an effect on the variances and 
covariance of components; the variance of one component is proportional to the force constant of the other and to 
the coupling strength between the two, and the covariance between the two is proportional to the coupling strength. 

We apply the above results to study the regulation of an auto-regulatory network of a single gene, which describes 
the central dogma of gene expression, transcription and translation. The two variables, x\ and x 2 , refer to the 
concentrations of mRNA and protein, respectively. In this study, we use the most common noise-attenuating regulatory 



Here, k 



max 



mechanism, called negative feedback and described by Hill function i?i (x 2 ) = fc max / 1 + (a^/fcd)' 3 

the maximum transcription rate of mRNA, kd is the binding constant specifying the threshold protein concentration 
at which the transcription rate is half its maximum value, and (3 is the Hill coefficient. On the other hand, we 
set R2 (xi) = fc 2 xi, where fc 2 is the translation rate of protein. Then, a stable equilibrium, x* , is characterized by 
two conditions: <fi 1 (j) 2 — r\ (x 2 ) k 2 > and <p 1 + 4> 2 > 0. Subsequently, one can use Bendixson's criterion to further 
conlude that there is no any cycles, only one equilibriuum exists [Hj]. For the values of the parameters, we mainly 
follow those given in Refs. [a El- The half-lifcs of mRNA molecules and proteins are set as 2 minutes and 1 hour, 
respectively; this leads to <fi 1 = (In 2) /2 and <f> 2 — (In 2) /60 in the unit of (min) . The average size of a burst of 
proteins, b — k 2 /(j> x , is set as 10, this leads to k 2 = 5 (ln2). By using the fact that the protein concentration is about 
1200 when (3 — (no feedback), we set fc max = 3[l(|. To study the effect of the strength of negative feedback on the 
characteristics of the system, we vary the parameters /3 from 2 to 11, while the kd value is fixed as 800. 

The numerical results are shown in Fig. 1 (a) for the K\ and k 2 values and in Fig. 1 (b) for the k\ 2 and AG values 
as functions of the equilibrium concentration of protein x 2 for different values. As a consequence of increasing the (3 
value, the k\ and n 2 values also increase but the k\ 2 value decreases; the values are ranged between 0.305 < n\ < 0.354, 
-6.00 x 10" 4 < k 2 < 6.26 x 10~ 4 , and 2.74 x 10~ 4 < k\ 2 < 8.08 x 10" 4 , respectively. Note that the n 2 values are 
drastically smaller than the K\ values, reflecting the longer half-life of protein. Moreover, the AG values for different 
(3 given in Fig. 1(b) indicate that the system becomes more in focus when the (3 value increases. 

The fluctuations of the system near a stable point are analyzed by measuring the variance of Xi in terms of the Fano 
factors Vi, defined as Vi — £1 (ct^/x*), and the covariance of x\ and x 2 in terms of the correlation coefficient i?i 2) 

denned as i?i 2 = a\ 2I \j a \,\ a \ 2- The numerical results of v\, v%, and i?i 2 for different (3 are shown in Fig. 2. The 
plots indicate that a larger k 2 implies a smaller z/ 2 and a larger n\ 2 implies a larger i?i 2 . Furthermore, as indicated in 
the inset of Fig. 2, the v x values all are very closed to but less than one over the range of 2 < {3 < 11; the value firstly 
decreases from V\ = 0.9844 at j3 = 2, reaches the minimum v\ = 0.9827 at j3 = 4, and then increases to v\ = 0.9926 at 
(3 = 11. Because that the fluctuation of mRNA is caused by a process very close to Poissonian with V\ = 1, we have 
very small correlation coefficient ranged between 0.015 < i?i 2 < 0.101. We further compare the results with those 
obtained from the linear noise Fokker-Planck equation which, as shown explicitly in Ref . [ll | , describes the distribution 
of fluctuations (i) introduced via the setting, x, (t) = Xi (t) + fi -1 / 2 ^ (t), where the macroscopic values 5, (f) are 
determined by the rate equations of Eq. (fT]). The results thus obtained agree completely with those shown in Fig. 2. 

In conclusion, we present a mechanical viewpoint on the characteristics of regulatory genetic networks, which is 
obtained from the energy landscape of the network near a stable point. The new approach is shown to be consistent 
with other descriptions, as demonstrated in the explicit calculations of auto-regulatory networks, it also provides 
additional informations, such as the steepness of a stable point and its relation to fluctuations. Though the method 
can also be applied to other genetic networks straightforwardly, it is limited in the sense that the overall potential 
landscape cannot be approximated by a local quadratic approach. For bi-stable systems such as toggle switches, one 
might need to patch the potential derived from the two local minimums in a rigorous fashion to obtain a better result. 
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FIG. 1: The force constants, ki, K2, and /t 1>2 , and the effective free energy difference AG for different /3 values: (a) ki (solid 
squares) and K2 (hollow squares), (b) (solid squares) and AG (hollow squares). The horizontal axis is the equilibrium 
concentration of protein x\ ■ 



FIG. 2: The Fano factors, v\ and h>2, and the correlation coefficient R12 for different (5 values. The horizontal axis is the 
equilibrium concentration of protein x?, and the inset shows the details of the v\ values. 



